//----------------------------------------------------------------------------
//
// (C) Copyrights MARIN/MSCN  2002 
//
//----------------------------------------------------------------------------
//
// $Author: A. Abbing $
//
// $Description: $
//
// $End: $
//
//----------------------------------------------------------------------------

#ifndef __mst_Vector3_h__
#define __mst_Vector3_h__

#include "Math.h"
#include "string.h"
#include "Constants.h"
#include "Matrix4.h"
#include "Stream.h"
#include <iostream>

#if defined(min) && !defined(mst_min)
#define mst_min min
#undef min
#endif

#if defined(max) && !defined(mst_max)
#define mst_max max
#undef max
#endif


namespace mst {

/*! 
  \class Vector3 mst/Vector3.h
  \brief 3D vector template class
  \ingroup math
*/

  template<class Type> class Vector3 {
      Type m_v[3];

    public:
      /** return the size of the vector */
      virtual unsigned int size() const { return 3; }

      /** Get vector internal data. */
      Type const* getData() const {
        return m_v;
      }

      /** Get vector internal data. */
      Type* getData() {
        return const_cast<Type*>( static_cast<Vector3 const&>(*this).getData() );
      }

      /*!
	\fn Vector3()
	\brief Construct a 3D vector with initial values (0.0, 0.0, 0.0)
      */
      Vector3() { memset((void*) m_v, 0, 3*sizeof(Type)); }
	
      /*!
	\fn Vector3(Type x, Type y, Type z)
	\brief Construct a 3D vector
	\param x the x initialisation value
	\param y the x initialisation value
	\param z the x initialisation value
      */
      Vector3(Type x, Type y, Type z) { set(x, y, z); }

      /*!
	\fn Vector3(Type v)
	\brief Construct a 3D vector
	\param v the x,y,z initialisation value
      */
      Vector3(Type v) { set(v, v, v); }

      /**Copy constructor */
      Vector3(const Vector3<Type> &v) { set(v[0], v[1], v[2]); }

      /*!
	\fn virtual ~Vector3()
	\brief Destruct a mst::Vector3 object.
      */
      virtual ~Vector3() {}

      /*!
	\fn void set(const Vector3<Type> &v)
	\param v the initialisation vector
	\brief Set the value of the vector
      */
      void set(const Vector3<Type> &v) { m_v[0] = v[0]; m_v[1] = v[1]; m_v[2] = v[2]; }
      /*!
	\fn void get(Type &x, Type &y, Type &z) const
	\brief Get the values of the vector
	\param x the x value
	\param y the y value
	\param z the z value
      */
      void get(Type &x, Type &y, Type &z) const 
      { x = m_v[0]; y = m_v[1]; z = m_v[2]; }
      /*!
	\fn void set(Type x, Type y, Type z)
	\brief Set the value of the vector
	\param x the x initialisation value
	\param y the y initialisation value
	\param z the z initialisation value
      */
      void set(Type x, Type y, Type z) { m_v[0] = x; m_v[1] = y; m_v[2] = z; }
	
      /*!
	\fn void clip(const Vector3<Type> &minv, const Vector3<Type> &maxv)
	\brief clip the vector against the minimum and maximum vector values
	\param minv the minimum clip value
	\param maxv the maximum clip value
      */
      void clip(const Vector3<Type> &minv, const Vector3<Type> &maxv) {
        m_v[0] = mst::clamp(m_v[0], minv[0], maxv[0]);
        m_v[1] = mst::clamp(m_v[1], minv[1], maxv[1]);
        m_v[2] = mst::clamp(m_v[2], minv[2], maxv[2]);
      }
	 
      /*!
	\fn Type sqrLength() const
	\brief return the squared length of the vector
      */
      Type sqrLength() const { 
        return m_v[0]*m_v[0] + m_v[1]*m_v[1] + m_v[2]*m_v[2]; 
      }
	
      /*!
	\fn Type length() const
	\brief return the length of the vector
      */
      Type length() const { 
        return mst::sqrt(sqrLength());
      }

      /*!
	\fn Type sqrDistance(const Vector3<Type> &v) const
	\param v the second vector
	\brief return the square distance between two vectors
      */
      Type sqrDistance(const Vector3<Type> &v) const {
        Type x = v[0] - m_v[0];
        Type y = v[1] - m_v[1];
        Type z = v[2] - m_v[2];
        return x*x + y*y + z*z;
      }
	
      /*!
	\fn Type distance(const Vector3<Type> &v) const
	\param v the second vector
	\brief return the distance between two vectors
      */
      Type distance(const Vector3<Type> &v) const { return mst::sqrt(sqrDistance(v)); }

      /*!
	\fn Type &operator[](int i)
	\param i index
	\brief Accesses indexed component of vector
      */
      Type &operator[](int i) { return m_v[i]; }

      /*!
	\fn const Type &operator[](int i) const
	\param i index
	\brief Accesses indexed component of vector
      */
      const Type &operator[](int i) const { return m_v[i]; }

      /*! \fn void max(const Vector3<Type> &m)
        \brief Maximize the this vector
        \param m The maximum value
      */
      void max(const Vector3<Type> &m)
      {
        m_v[0] = mst::max(m_v[0], m[0]);
        m_v[1] = mst::max(m_v[1], m[1]);
        m_v[2] = mst::max(m_v[2], m[2]);
      }

      /*! \fn void max(const Vector3<Type> &v, const Vector3<Type> &m)
        \brief Maximize the this vector 
        \param v The value vector
        \param m The maximum value
      */
      void max(const Vector3<Type> &v, const Vector3<Type> &m)
      {
        m_v[0] = mst::max(v[0], m[0]);
        m_v[1] = mst::max(v[1], m[1]);
        m_v[2] = mst::max(v[2], m[2]);
      }

      /*! \fn void min(const Vector3<Type> &m)
        \brief Minimize the this vector
        \param m The minimum value
      */
      void min(const Vector3<Type> &m)
      {
        m_v[0] = mst::min(m_v[0], m[0]);
        m_v[1] = mst::min(m_v[1], m[1]);
        m_v[2] = mst::min(m_v[2], m[2]);
      }

      /*! \fn void min(const Vector3<Type> &v, const Vector3<Type> &m)
        \brief Minimize the this vector 
        \param v The value vector
        \param m The minimum value
      */
      void min(const Vector3<Type> &v, const Vector3<Type> &m)
      {
        m_v[0] = mst::min(v[0], m[0]);
        m_v[1] = mst::min(v[1], m[1]);
        m_v[2] = mst::min(v[2], m[2]);
      }

      /*!
	\fn void copy(const Vector3<Type> &v)
	\param v Vector to copy
	\brief Copy vector
      */
	
      void copy(const Vector3<Type> &v) { set(v[0], v[1], v[2]); }
      /*!
	\fn void negate(const Vector3<Type> &v)
	\brief Negate the current vector
      */
	
      void negate(const Vector3<Type> &v) { set(-v[0], -v[1], -v[2]); }
      /*!
	\fn Type normalize()
	\brief Normalise the current vector
      */
	
      Type normalize() {
        Type l = length();
        *this *= (l < mst::eta(l)) ? mst::huge(l) : ((Type)1) / l;
        return l;
      }
	
      /*!
	\fn void scale(Type s, const Vector3<Type> &v)
	\param s The fraction of the vector
	\param v The vector
	\brief Set a scaled vector
      */
      void scale(Type s, const Vector3<Type> &v) { set(s*v[0], s*v[1], s*v[2]); }
	
      /*!
	\fn void addScaled(const Vector3<Type>& a, Type s, const Vector3<Type>& b)
	\param a The first vector
	\param s The fraction of vector b
	\param b The second vector
	\brief Add a scaled vector
      */
      void addScaled(const Vector3<Type>& a, Type s, const Vector3<Type>& b) {
        m_v[0] = a[0] + s * b[0];
        m_v[1] = a[1] + s * b[1];
        m_v[2] = a[2] + s * b[2];
      }
	
      /*!
	\fn void subScaled(const Vector3<Type>& a, Type s, const Vector3<Type>& b)
	\param a The first vector
	\param s The fraction of vector b
	\param b The second vector
	\brief Subtract a scaled vector
      */
      void subScaled(const Vector3<Type>& a, Type s, const Vector3<Type>& b) {
        m_v[0] = a[0] - s * b[0];
        m_v[1] = a[1] - s * b[1];
        m_v[2] = a[2] - s * b[2];
      }
	
      /*!
	\fn void combine(Type sa, const Vector3<Type>& a, Type sb, const Vector3<Type>& b)
	\param sa The fraction of vector a
	\param a The first vector
	\param sb The fraction of vector b
	\param b The second vector
	\brief Combine two vectors
      */
      void combine(Type sa, const Vector3<Type>& a, Type sb, const Vector3<Type>& b) {
        m_v[0] = sa * a[0] + sb * b[0];
        m_v[1] = sa * a[1] + sb * b[1];
        m_v[2] = sa * a[2] + sb * b[2];
      }
      /*!
	\fn bool equal(const Vector3<Type> &v) const
	\param v The vector
	\brief Test for component-wise equality
      */
      bool equal(const Vector3<Type> &v) const {
        return m_v[0] == v[0] && m_v[1] == v[1] && m_v[2] == v[2];
      }
	
      /*!
	\fn bool almostEqual(const Vector3<Type>& v, Type tol) const
	\param v the second vector
	\param tol the tolerance value
	\brief Tests for approximate component-wise equality
	\sa equal

	Tests for approximate component-wise equality of
	two vectors. It returns false or true depending on whether
	the absolute value of the difference between each pair of components is
	less than the tolerance. 
      */
      bool almostEqual(const Vector3<Type>& v, Type tol) const {
        return mst::equal(m_v[0], v[0], tol) &&
          mst::equal(m_v[1], v[1], tol) &&
          mst::equal(m_v[2], v[2], tol);
      }

      /*!
	\fn bool equal(const Vector3<Type>& v, Type tol) const
	\param v the second vector
	\param tol the tolerance value
	\brief Tests for approximate component-wise equality
	\sa almostEqual

	Tests for approximate component-wise equality of
	two vectors. It returns false or true depending on whether
	the absolute value of the difference between each pair of components is
	less than the tolerance. 
      */
      bool equal(const Vector3<Type>& v, Type tol) const 
      { return almostEqual(v, tol); }

      /*!
	\fn bool almostEqual(const Vector3<Type>& v, const Vector3<Type>& tol) const
	\param v the second vector
	\param tol the tolerance vector
	\brief Tests for approximate component-wise equality
	\sa equal

	Tests for approximate component-wise equality of
	two vectors. It returns false or true depending on whether
	the absolute value of the difference between each pair of components is
	less than the tolerance. 
      */
      bool almostEqual(const Vector3<Type>& v, const Vector3<Type>& tol) const {
        return mst::equal(m_v[0], v[0], tol[0]) &&
          mst::equal(m_v[1], v[1], tol[1]) &&
          mst::equal(m_v[2], v[2], tol[2]);
      }
	
      /*!
	\fn bool equal(const Vector3<Type>& v, const Vector3<Type>& tol) const
	\param v the second vector
	\param tol the tolerance vector
	\brief Tests for approximate component-wise equality
	\sa almostEqual

	Tests for approximate component-wise equality of
	two vectors. It returns false or true depending on whether
	the absolute value of the difference between each pair of components is
	less than the tolerance. 
      */
      bool equal(const Vector3<Type>& v, const Vector3<Type>& tol) const 
      { return almostEqual(v, tol); }
	
      /*!
	\fn void cross(const Vector3<Type>& a, const Vector3<Type>& b)
	\param a the first vector
	\param b the second vector
	\brief Calculate the cross product of two vectors
      */
      void cross(const Vector3<Type>& a, const Vector3<Type>& b) {
        m_v[0] = a[1]*b[2] - a[2]*b[1];
        m_v[1] = a[2]*b[0] - a[0]*b[2];
        m_v[2] = a[0]*b[1] - a[1]*b[0];
      }
	
      /*!
	\fn Type dot(const Vector3<Type>& v) const
	\param v the second vector
	\brief Calculate the dot product of two vectors
      */
      Type dot(const Vector3<Type>& v) const {
        return m_v[0]*v[0] + m_v[1]*v[1] + m_v[2]*v[2];
      }

      /*!
	\fn void add(const Vector3<Type> &a, const Vector3<Type> &b)
	\param a The first vector
	\param b The second vector
	\brief Add two vectors
      */
      void add(const Vector3<Type> &a, const Vector3<Type> &b) {
        m_v[0] = a[0] + b[0];
        m_v[1] = a[1] + b[1];
        m_v[2] = a[2] + b[2];
      }

      /*!
	\fn void sub(const Vector3<Type> &a, const Vector3<Type> &b)
	\param a The first vector
	\param b The second vector
	\brief Subtract two vectors
      */
      void sub(const Vector3<Type> &a, const Vector3<Type> &b) {
        m_v[0] = a[0] - b[0];
        m_v[1] = a[1] - b[1];
        m_v[2] = a[2] - b[2];
      }

      /*!
	\fn void xformPt(const Vector3<Type>& v, const mst::Matrix4<Type>& m)
	\param v The vector
	\param m The matrix
	\brief Transform of a point
      */
      void xformPt(const Vector3<Type>& v, const mst::Matrix4<Type>& m) {
        Type t0, t1, t2, w;
        t0 = v[0];
        t1 = v[1];
        t2 = v[2];
        m_v[0] = (t0 * m[0][0] + t1 * m[1][0] + t2 * m[2][0] + m[3][0]);
        m_v[1] = (t0 * m[0][1] + t1 * m[1][1] + t2 * m[2][1] + m[3][1]);
        m_v[2] = (t0 * m[0][2] + t1 * m[1][2] + t2 * m[2][2] + m[3][2]);
        w = (t0*m[0][3] + t1*m[1][3] + t2*m[2][3] + m[3][3]);
        if (w != (Type) 1) {
          if (mst::abslt(w, mst::eta(w))) w = mst::eta(w);
          *this /= w;
        }
      }

      /*!
	\fn void xformPt(const mst::Matrix4<Type>& m)
	\param m The matrix
	\brief Transform of a point
      */
      void xformPt(const mst::Matrix4<Type>& m) {
        Type t0, t1, t2, w;
        t0 = m_v[0];
        t1 = m_v[1];
        t2 = m_v[2];
        m_v[0] = (t0 * m[0][0] + t1 * m[1][0] + t2 * m[2][0] + m[3][0]);
        m_v[1] = (t0 * m[0][1] + t1 * m[1][1] + t2 * m[2][1] + m[3][1]);
        m_v[2] = (t0 * m[0][2] + t1 * m[1][2] + t2 * m[2][2] + m[3][2]);
        w = (t0*m[0][3] + t1*m[1][3] + t2*m[2][3] + m[3][3]);
        if (w != (Type) 1) {
          if (mst::abslt(w, mst::eta(w))) w = mst::eta(w);
          *this /= w;
        }
      }

      /*!
	\fn void xformVec(const Vector3<Type>& v, const mst::Matrix4<Type>& m)
	\param v The vector
	\param m The matrix
	\brief Transform of a vector
      */
      void xformVec(const Vector3<Type>& v, const mst::Matrix4<Type>& m) {
        Type t0, t1, t2;
        t0 = v[0]; t1 = v[1]; t2 = v[2];
        m_v[0] = (t0 * m[0][0] + t1 * m[1][0] + t2 * m[2][0]);
        m_v[1] = (t0 * m[0][1] + t1 * m[1][1] + t2 * m[2][1]);
        m_v[2] = (t0 * m[0][2] + t1 * m[1][2] + t2 * m[2][2]);
      }

      /*!
	\fn void xformVec(const mst::Matrix4<Type>& m)
	\param m The matrix
	\brief Transform of a vector
      */
      void xformVec(const mst::Matrix4<Type>& m) {
        Type t0, t1, t2;
        t0 = m_v[0]; t1 = m_v[1]; t2 = m_v[2];
        m_v[0] = (t0 * m[0][0] + t1 * m[1][0] + t2 * m[2][0]);
        m_v[1] = (t0 * m[0][1] + t1 * m[1][1] + t2 * m[2][1]);
        m_v[2] = (t0 * m[0][2] + t1 * m[1][2] + t2 * m[2][2]);
      }

      /*!
	\fn void fullXformPt(const Vector3<Type>& v, const mst::Matrix4<Type>& m)
	\param v The vector
	\param m The matrix
	\brief Full transform of a vector
      */
      void fullXformPt(const Vector3<Type>& v, const mst::Matrix4<Type>& m){
        Type t0, t1, t2, w;
        t0 = v[0];
        t1 = v[1];
        t2 = v[2];
        m_v[0] = (t0 * m[0][0] + t1 * m[1][0] + t2 * m[2][0] + m[3][0]);
        m_v[1] = (t0 * m[0][1] + t1 * m[1][1] + t2 * m[2][1] + m[3][1]);
        m_v[2] = (t0 * m[0][2] + t1 * m[1][2] + t2 * m[2][2] + m[3][2]);
        w = (t0*m[0][3] + t1*m[1][3] + t2*m[2][3] + m[3][3]);
        if (w != (Type) 1) {
          if (mst::abslt(w, mst::eta(w))) w = mst::eta(w);
          *this /= w;
        }
      }

      /*!
	\fn void fullXformPt(const mst::Matrix4<Type>& m)
	\param m The matrix
	\brief Full transform of a vector
      */
      void fullXformPt(const mst::Matrix4<Type>& m){
        Type t0, t1, t2, w;
        t0 = m_v[0];
        t1 = m_v[1];
        t2 = m_v[2];
        m_v[0] = (t0 * m[0][0] + t1 * m[1][0] + t2 * m[2][0] + m[3][0]);
        m_v[1] = (t0 * m[0][1] + t1 * m[1][1] + t2 * m[2][1] + m[3][1]);
        m_v[2] = (t0 * m[0][2] + t1 * m[1][2] + t2 * m[2][2] + m[3][2]);
        w = (t0*m[0][3] + t1*m[1][3] + t2*m[2][3] + m[3][3]);
        if (w != (Type) 1) {
          if (mst::abslt(w, mst::eta(w))) w = mst::eta(w);
          *this /= w;
        }
      }

      /*!
	\fn void angles(const Vector3<Type>& v, Type &pitch, Type &yaw) const
	\param v The target vector
	\param pitch the pitch angle of the vector
	\param yaw the yaw angle of the vector
	\brief Calculate the pitch and yaw angle of the vector
      */
      void angles(const Vector3<Type>& v, Type &pitch, Type &yaw) const {
        Vector3<Type> a(m_v[0]-v[0], m_v[1]-v[1], m_v[2]-v[2]);
        pitch = mst::atan2(a[2], mst::sqrt(a[0]*a[0] + a[1]*a[1]));
        yaw = mst::atan2(a[1], a[0]);
      }

      /*!
	\fn void angles(Type &pitch, Type &yaw) const
	\param pitch the pitch angle of the vector
	\param yaw the yaw angle of the vector
	\brief Calculate the pitch and yaw angle of the vector
      */
      void angles(Type &pitch, Type &yaw) const {
        pitch = mst::atan2(m_v[2], mst::sqrt(m_v[0]*m_v[0] + m_v[1]*m_v[1]));
        yaw = mst::atan2(m_v[1], m_v[0]);
      }

      /*!
	\fn bool operator==(const Vector3<Type> &v)
	\param v The vector
	\brief Test for equality of a vector with the current vector
      */
      bool operator==(const Vector3<Type> &v) { return equal(v); }

      /*!
	\fn bool operator!=(const Vector3<Type> &v)
	\param v The vector
	\brief Test for inequality of a vector with the current vector
      */
      bool operator!=(const Vector3<Type> &v) { return !equal(v); }

      /*!
	\fn void operator-=(const Vector3<Type> &v)
	\param v The vector
	\brief Subtract a vector of the current vector
      */
      void operator-=(const Vector3<Type> &v) {
        m_v[0] -= v[0];
        m_v[1] -= v[1];
        m_v[2] -= v[2];
      }

      /*!
	\fn void operator+=(const Vector3<Type> &v)
	\param v The vector
	\brief Add a vector to the current vector
      */
      void operator+=(const Vector3<Type> &v) {
        m_v[0] += v[0];
        m_v[1] += v[1];
        m_v[2] += v[2];
      }

      /*!
	\fn void operator/=(Type s)
	\param s The scalar
	\brief Devide the current vector by a scalar
      */
      void operator/=(Type s) {
        m_v[0] /= s;
        m_v[1] /= s;
        m_v[2] /= s;
      }

      /*!
	\fn void operator*=(Type s)
	\param s The scalar
	\brief Multiply the current vector by a scalar
      */
      void operator*=(Type s) {
        m_v[0] *= s;
        m_v[1] *= s;
        m_v[2] *= s;
      }
      /*!
	\fn void operator*=(const Matrix4<Type> &m)
	\param m The 4x4 matrix
	\brief Multiply the current vector by a matrix
      */
      void operator*=(const Matrix4<Type> &m)
      {
        // transform as point (w=1), assuming affine transformation
        // i.e. does not use slower dst.xformFullPt().
        this->xformPt(*this, m);
      }

      /*! 
	\fn bool operator<(const Vector3<Type> &v)
	\brief Componentwise smaller than check
      */
      bool operator<(const Vector3<Type> &v)
      {
        return (m_v[0] < v[0]) || (m_v[1] < v[1]) || (m_v[2] < v[2]);
      }

      /*! 
	\fn bool operator>(const Vector3<Type> &v)
	\brief Componentwise larger than check
      */
      bool operator>(const Vector3<Type> &v)
      {
        return (m_v[0] > v[0]) || (m_v[1] > v[1]) || (m_v[2] > v[2]);
      }

    public:	// Warning: return by value can be slow
      /*!
	\fn Vector3<Type> operator-() const
	\brief Return the result of a unary vector negation
      */
      Vector3<Type> operator-() const { return Vector3<Type>(-m_v[0], -m_v[1], -m_v[2]); }
      /*!
	\fn Vector3<Type> operator+(const Vector3<Type>& v) const
	\param v The 3D vector
	\brief Return the result of  a vector addition
      */
      Vector3<Type> operator+(const Vector3<Type>& v) const {
        return Vector3<Type>(m_v[0]+v[0], m_v[1]+v[1], m_v[2]+v[2]);
      }
      /*!
	\fn Vector3<Type> operator-(const Vector3<Type>& v) const
	\param v The 3D vector
	\brief Return the result of a vector subtraction
      */
      Vector3<Type> operator-(const Vector3<Type>& v) const {
        return Vector3<Type>(m_v[0]-v[0], m_v[1]-v[1], m_v[2]-v[2]);
      }
      /*!
	\fn Vector3<Type>& operator=(const Vector3<Type>& v)
	\param v The 3D vector
	\brief Assign a vector 
      */
      Vector3<Type>& operator=(const Vector3<Type>& v) {
        set(v[0], v[1], v[2]);
        return *this;
      }
      /*!
	\fn friend Vector3<Type> operator*(Type s, const Vector3<Type>& v)
	\param s The scalar
	\param v The 3D vector
	\brief Multiply a vector by a scalar
      */
      friend Vector3<Type> operator*(Type s, const Vector3<Type>& v)
      { return Vector3<Type>(v[0]*s, v[1]*s, v[2]*s); }
      /*!
	\fn friend Vector3<Type> operator*(const Vector3<Type>& v, Type s)
	\param v The 3D vector
	\param s The scalar
	\brief Multiply a vector by a scalar
      */
      friend Vector3<Type> operator*(const Vector3<Type>& v, Type s)
      { return Vector3<Type>(v[0]*s, v[1]*s, v[2]*s); }
      /*!
	\fn friend Vector3<Type> operator/(const Vector3<Type>& v, Type s)
	\param v The 3D vector
	\param s The scalar
	\brief Devide a vector by a scalar
      */
      friend Vector3<Type> operator/(const Vector3<Type>& v, Type s)
      { 	s = ((Type)1)/s; 
      return Vector3<Type>(v[0]*s, v[1]*s, v[2]*s);
      }
      /*!
	\fn friend Vector3<Type> operator*(const Vector3<Type>& v, const mst::Matrix4<Type>& m)
	\param v The 3D vector
	\param m The 4x4 homogeneous matrix
	\brief multiply a vector by a matrix
      */
      friend Vector3<Type> operator*(const Vector3<Type>& v, const mst::Matrix4<Type>& m)
      {
        // transform as point (w=1), assuming affine transformation
        // i.e. does not use slower dst.xformFullPt().
        Vector3<Type> dst; dst.xformPt(v, m); return dst;
      }

      /*!
	\fn friend Vector3<Type> operator*(const Vector3<Type>& a, const Vector3<Type>& b)
	\param a The 3D vector
	\param b The 3D vector
	\brief multiply 2 3D vectors
      */
      friend Vector3<Type> operator*(const Vector3<Type>& a, const Vector3<Type>& b)
      {	Vector3<Type> r(a[0]*b[0], a[1]*b[1], a[2]*b[2]);
      return r;
      }
  };

/*! 
  \typedef Vector3f
  \brief float specialisation of a 3D vector
  \ingroup math
*/

  typedef Vector3<float> Vector3f;

/*! 
  \typedef Vector3d
  \brief double specialisation of a 3D vector
  \ingroup math
*/
  typedef Vector3<double> Vector3d;

/*! 
  \fn template<class Type> std::ostream& operator<<(std::ostream &s, const mst::Vector3<Type> &v)
  \brief stream output operator of a 3D vector
  \ingroup io
*/
  template<class Type> std::ostream& operator<<(std::ostream &s, const mst::Vector3<Type> &v)
  { s << v[0] << " " << v[1] << " " << v[2]; return s; }

/*! 
  \fn template<class Type> std::istream& operator>>(std::istream &s, mst::Vector3<Type> &v)
  \brief stream input operator of a 3D vector
  \ingroup io
*/
  template<class Type> std::istream& operator>>(std::istream &s, mst::Vector3<Type> &v)
  { char c; s >> v[0] >> c >> v[1] >> c >> v[2]; return s; }

/*! 
  \fn template<class T> mst::Vector3<T> abs(const mst::Vector3<T> &v)
  \brief return vector with absolute values
  \ingroup math
*/
  template<class T> mst::Vector3<T> abs(const mst::Vector3<T> &v)
  {	mst::Vector3<T> res(mst::abs(v[0]), mst::abs(v[1]), mst::abs(v[2]));
  return res;
  }

/*! 
  \fn template<class Type> mst::ostream& operator<<(mst::ostream &s, const mst::Vector3<Type> &v)
  \brief mst::ostream output operator
  \ingroup io
*/
  template<class Type> mst::ostream& operator<<(mst::ostream &s,
                                                const mst::Vector3<Type> &v)
  {
    s << v[0] << v[1] << v[2];
    return s;
  }

/*! 
  \fn template<class Type> mst::istream& operator>>(mst::istream &s, mst::Vector3<Type> &v)
  \brief mst::istream input operator
  \ingroup io
*/
  template<class Type> mst::istream& operator>>(mst::istream &s, 
                                                mst::Vector3<Type> &v)
  {
    s >> v[0] >> v[1] >> v[2];
    return s;
  }


}

#endif
